#library(soilP)
library(here)
library(devtools)
load_all("./")
Loading bugcount
data_dir <- file.path(system.file(package = "bugcount"),"extdata")

1 Whitefly analysis.

1.1 Read WF count data.

WF counts per picture. Each picture must have a unique combination of
experiment + replicate + plant + pic
experiment description consists of:
exp_year + exp_cross + exp_propagation + exp_substrate

wf_file <- file.path(data_dir, "wf_consolidated.tab")
wf <- read_wf(file = wf_file, sep = "\t",header=TRUE)

1.2 Counts per Plant

wf_count <- plant_count(wf)
summary(wf_count$nymphs)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0     555    1859    2506    3460   25317 

What indicator should we use as resistance phenotype? Nymph counts are not normal Error is not distributed normally. Nymph counts do not meet linear model assumptions (ANOVA, MapQTL). Best fit is a negative Binomial Distribution

1.3 Probability density model fit to actual distribution

wf_count_fit <-  count_fit(wf_count)
plot_count_fit(wf_count_fit, main = "Nymphs")

1.4 Reproducibility

Assume that nymphs per plant are distributed as negative binomial

1.4.1 Whole population.

nymphs ~ experiment GLM, posthoc and common letter difference

by_experiment <- fit_nb_glm(nymphs ~ experiment, wf_count)
plot_fit_nb_glm(wf_count,by_experiment, type = "violin", cld = FALSE)

add infestation regime as factor high,low: inferred from nyphs ~ experiment GLM

wf_count$infestation <- "low"
wf_count[wf_count$exp_year > 2016,]$infestation <- "high"
wf_count$infestation <- factor(wf_count$infestation, levels =c("low","high"))

1.5 Correlation between experiments

cor_title <- paste("Geometric Means Correlation","Log Scale", sep = "\n")
for (cross in levels(wf_count$exp_cross)) {
  
  wf_by_exp <- wf_wide(clone ~ experiment,
                       wf_count[wf_count$exp_cross == cross,])
  
  plot_wide_cor(wf_by_exp, main = paste(cross,cor_title))
  
}

1.6 Checks.

wf_check <- wf_count[grepl("internal_check", wf_count$group) &
                       wf_count$clone != "Secundina",]
wf_check <- remove_levels(wf_check)
checks_count_fit <-  count_fit(wf_check)
plot_count_fit(checks_count_fit, main = "Nymphs")

1.6.1 Correlation between experiments

check_by_exp <- wf_wide(clone ~ experiment, wf_check, 
                        fun = function(x) log10(geometric_mean(x)) )
plot_check_cor(check_by_exp, main = paste("Checks", cor_title))

# nymphs ~ clone GLM posthoc and common letter difference
fit <- fit_nb_glm(nymphs ~ clone, wf_check)
plot_fit_nb_glm(wf_check,fit,type = "violin", xmax = 10000)

2 Nymph counts by Cross

LS0tCnRpdGxlOiAiQnVnIENvdW50IEFuYWx5c2lzIgpvdXRwdXQ6CiAgaHRtbF9ub3RlYm9vazoKICAgIGhpZ2hsaWdodDogdGFuZ28KICAgIG51bWJlcl9zZWN0aW9uczogeWVzCiAgICB0aGVtZTogc3BhY2VsYWIKICAgIHRvYzogeWVzCiAgICB0b2NfZmxvYXQ6IHllcwotLS0KCmBgYHtyfQojbGlicmFyeShzb2lsUCkKbGlicmFyeShoZXJlKQpsaWJyYXJ5KGRldnRvb2xzKQpsb2FkX2FsbCgiLi8iKQpkYXRhX2RpciA8LSBmaWxlLnBhdGgoc3lzdGVtLmZpbGUocGFja2FnZSA9ICJidWdjb3VudCIpLCJleHRkYXRhIikKCmBgYAoKIyBXaGl0ZWZseSBhbmFseXNpcy4KIyMgUmVhZCBXRiBjb3VudCBkYXRhLiAKV0YgY291bnRzIHBlciBwaWN0dXJlLiBFYWNoIHBpY3R1cmUgbXVzdCBoYXZlIGEgdW5pcXVlIGNvbWJpbmF0aW9uIG9mICAKZXhwZXJpbWVudCArIHJlcGxpY2F0ZSArIHBsYW50ICsgcGljICAKZXhwZXJpbWVudCBkZXNjcmlwdGlvbiBjb25zaXN0cyBvZjogIApleHBfeWVhciArIGV4cF9jcm9zcyArIGV4cF9wcm9wYWdhdGlvbiArIGV4cF9zdWJzdHJhdGUgIAoKCmBgYHtyfQp3Zl9maWxlIDwtIGZpbGUucGF0aChkYXRhX2RpciwgIndmX2NvbnNvbGlkYXRlZC50YWIiKQp3ZiA8LSByZWFkX3dmKGZpbGUgPSB3Zl9maWxlLCBzZXAgPSAiXHQiLGhlYWRlcj1UUlVFKQpgYGAKIyMgQ291bnRzIHBlciBQbGFudCAgCgpgYGB7cn0KCndmX2NvdW50IDwtIHBsYW50X2NvdW50KHdmKQpzdW1tYXJ5KHdmX2NvdW50JG55bXBocykKCmBgYAoKV2hhdCBpbmRpY2F0b3Igc2hvdWxkIHdlIHVzZSBhcyByZXNpc3RhbmNlIHBoZW5vdHlwZT8KTnltcGggY291bnRzIGFyZSBub3Qgbm9ybWFsIApFcnJvciBpcyBub3QgZGlzdHJpYnV0ZWQgbm9ybWFsbHkuCk55bXBoIGNvdW50cyBkbyBub3QgbWVldCBsaW5lYXIgbW9kZWwgYXNzdW1wdGlvbnMgKEFOT1ZBLCBNYXBRVEwpLgpCZXN0IGZpdCBpcyBhIG5lZ2F0aXZlIEJpbm9taWFsIERpc3RyaWJ1dGlvbgoKIyMgUHJvYmFiaWxpdHkgZGVuc2l0eSBtb2RlbCBmaXQgdG8gYWN0dWFsIGRpc3RyaWJ1dGlvbgpgYGB7ciwgd2FybmluZz1GQUxTRX0Kd2ZfY291bnRfZml0IDwtICBjb3VudF9maXQod2ZfY291bnQpCgpwbG90X2NvdW50X2ZpdCh3Zl9jb3VudF9maXQsIG1haW4gPSAiTnltcGhzIikKYGBgCgojIyBSZXByb2R1Y2liaWxpdHkgCgpBc3N1bWUgdGhhdCBueW1waHMgcGVyIHBsYW50IGFyZSBkaXN0cmlidXRlZCBhcyBuZWdhdGl2ZSBiaW5vbWlhbAoKIyMjIFdob2xlIHBvcHVsYXRpb24uCm55bXBocyB+IGV4cGVyaW1lbnQgR0xNLCBwb3N0aG9jIGFuZCBjb21tb24gbGV0dGVyIGRpZmZlcmVuY2UgCgpgYGB7cn0KYnlfZXhwZXJpbWVudCA8LSBmaXRfbmJfZ2xtKG55bXBocyB+IGV4cGVyaW1lbnQsIHdmX2NvdW50KQoKYGBgCgpgYGB7cn0KcGxvdF9maXRfbmJfZ2xtKHdmX2NvdW50LGJ5X2V4cGVyaW1lbnQsIHR5cGUgPSAidmlvbGluIiwgY2xkID0gRkFMU0UpCmBgYAphZGQgaW5mZXN0YXRpb24gcmVnaW1lIGFzIGZhY3RvciAKaGlnaCxsb3c6IGluZmVycmVkIGZyb20gbnlwaHMgfiBleHBlcmltZW50IEdMTQoKYGBge3J9CndmX2NvdW50JGluZmVzdGF0aW9uIDwtICJsb3ciCndmX2NvdW50W3dmX2NvdW50JGV4cF95ZWFyID4gMjAxNixdJGluZmVzdGF0aW9uIDwtICJoaWdoIgp3Zl9jb3VudCRpbmZlc3RhdGlvbiA8LSBmYWN0b3Iod2ZfY291bnQkaW5mZXN0YXRpb24sIGxldmVscyA9YygibG93IiwiaGlnaCIpKQpgYGAKCiMjIENvcnJlbGF0aW9uICBiZXR3ZWVuIGV4cGVyaW1lbnRzIAoKYGBge3IsIGZpZy5hc3A9MX0KCmNvcl90aXRsZSA8LSBwYXN0ZSgiR2VvbWV0cmljIE1lYW5zIENvcnJlbGF0aW9uIiwiTG9nIFNjYWxlIiwgc2VwID0gIlxuIikKCmZvciAoY3Jvc3MgaW4gbGV2ZWxzKHdmX2NvdW50JGV4cF9jcm9zcykpIHsKICAKICB3Zl9ieV9leHAgPC0gd2Zfd2lkZShjbG9uZSB+IGV4cGVyaW1lbnQsCiAgICAgICAgICAgICAgICAgICAgICAgd2ZfY291bnRbd2ZfY291bnQkZXhwX2Nyb3NzID09IGNyb3NzLF0pCiAgCiAgcGxvdF93aWRlX2Nvcih3Zl9ieV9leHAsIG1haW4gPSBwYXN0ZShjcm9zcyxjb3JfdGl0bGUpKQogIAp9CmBgYAoKIyMgQ2hlY2tzLgpgYGB7ciwgd2FybmluZz1GQUxTRX0Kd2ZfY2hlY2sgPC0gd2ZfY291bnRbZ3JlcGwoImludGVybmFsX2NoZWNrIiwgd2ZfY291bnQkZ3JvdXApICYKICAgICAgICAgICAgICAgICAgICAgICB3Zl9jb3VudCRjbG9uZSAhPSAiU2VjdW5kaW5hIixdCgp3Zl9jaGVjayA8LSByZW1vdmVfbGV2ZWxzKHdmX2NoZWNrKQoKY2hlY2tzX2NvdW50X2ZpdCA8LSAgY291bnRfZml0KHdmX2NoZWNrKQoKcGxvdF9jb3VudF9maXQoY2hlY2tzX2NvdW50X2ZpdCwgbWFpbiA9ICJOeW1waHMiKQpgYGAKCgojIyMgQ29ycmVsYXRpb24gIGJldHdlZW4gZXhwZXJpbWVudHMKCmBgYHtyLCBmaWcuYXNwPTF9CmNoZWNrX2J5X2V4cCA8LSB3Zl93aWRlKGNsb25lIH4gZXhwZXJpbWVudCwgd2ZfY2hlY2ssIAogICAgICAgICAgICAgICAgICAgICAgICBmdW4gPSBmdW5jdGlvbih4KSBsb2cxMChnZW9tZXRyaWNfbWVhbih4KSkgKQpwbG90X2NoZWNrX2NvcihjaGVja19ieV9leHAsIG1haW4gPSBwYXN0ZSgiQ2hlY2tzIiwgY29yX3RpdGxlKSkKYGBgCgoKYGBge3J9CiMgbnltcGhzIH4gY2xvbmUgR0xNIHBvc3Rob2MgYW5kIGNvbW1vbiBsZXR0ZXIgZGlmZmVyZW5jZQpmaXQgPC0gZml0X25iX2dsbShueW1waHMgfiBjbG9uZSwgd2ZfY2hlY2spCnBsb3RfZml0X25iX2dsbSh3Zl9jaGVjayxmaXQsdHlwZSA9ICJ2aW9saW4iLCB4bWF4ID0gMTAwMDApCmBgYAoKYGBge3IsIGZpZy5hc3A9MC41fQpmaXQgPC0gTUFTUzo6Z2xtLm5iKG55bXBocyB+ICBpbmZlc3RhdGlvbiAqIGNsb25lLCBkYXRhID0gd2ZfY2hlY2spCnBvc3Rob2MgPC0gbHNtZWFuczo6Y2xkKGxzbWVhbnM6OmxzbWVhbnMoZml0LCB+IGluZmVzdGF0aW9uICogY2xvbmUsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgYWRqdXN0ID0gInR1Y2tleSIpLCB0eXBlID0gInJlc3BvbnNlIikKIyBQbG90IHBvc3QgaG9jCgojIyMgT3JkZXIgdGhlIGxldmVscyBmb3IgcHJpbnRpbmcKY2xvbmVfb3JkZXIgPC0gd2ZfYWdncmVnYXRlKHdmX2NoZWNrLCBieV9zdGF0ID0gIm1lYW4iKVsiY2xvbmUiXQoKcG9zdGhvYyRjbG9uZSA8LSBmYWN0b3IocG9zdGhvYyRjbG9uZSwgbGV2ZWxzID0gY2xvbmVfb3JkZXIkY2xvbmUpCgpwb3N0aG9jJGluZmVzdGF0aW9uID0gZmFjdG9yKHBvc3Rob2MkaW5mZXN0YXRpb24sCiAgICAgICAgICAgICAgICAgICBsZXZlbHM9YygibG93IiwgImhpZ2giKSkKCiMjIyAgUmVtb3ZlIHNwYWNlcyBpbiAuZ3JvdXAgIAoKcG9zdGhvYyQuZ3JvdXA9Z3N1YigiICIsICIiLCBwb3N0aG9jJC5ncm91cCkKCiMgcXVhcnR6KCkKcGxvdF9nZ19pbmZlc3RhdGlvblhjbG9uZShwb3N0aG9jKQpgYGAKCgoKIyBOeW1waCBjb3VudHMgYnkgQ3Jvc3MgCmBgYHtyLCBmaWcuYXNwPTF9CmV4cF9sZXZlbHMgPC0gbGV2ZWxzKGFzLmZhY3Rvcih3Zl9jb3VudCRleHBlcmltZW50KSkKZXhwX2dycCA8LSBsaXN0KENNODk5NiA9IGxpc3QoMSwyLDQsNSksCiAgICAgICAgICAgICAgICBHTTg1ODYgPSBsaXN0KDMsNiw3KSkKcGxvdF9saXN0IDwtIGxpc3QoKQpmb3IgKGNyb3NzIGluIGxldmVscyh3Zl9jb3VudCRleHBfY3Jvc3MpKSB7CiAgZm9yIChpZHggaW4gMTpsZW5ndGgoZXhwX2dycFtbY3Jvc3NdXSkpIHsKICBleHBfYWxsb3dlZCA8LSBleHBfbGV2ZWxzW2V4cF9ncnBbW2Nyb3NzXV1bW2lkeF1dXQogIAogIHdmX2J5X2Nyb3NzIDwtIHdmX2NvdW50W3dmX2NvdW50JGV4cF9jcm9zcyA9PSBjcm9zcyAmCiAgICAgICAgICAgICAgICAgICAgICAgICAgICB3Zl9jb3VudCRleHBlcmltZW50ICVpbiUgZXhwX2FsbG93ZWQsXQogIGV4cF9jb3VudCA8LSBhZ2dyZWdhdGUoZXhwZXJpbWVudCB+IGNsb25lICx3Zl9ieV9jcm9zcywKICAgICAgICAgICAgICAgICAgICAgICAgIEZVTiA9IGZ1bmN0aW9uKHgpIGxlbmd0aCh1bmlxdWUoeCkpKQogIHNlbGVjdGVkX2Nsb25lcyA8LSB3aXRoKGV4cF9jb3VudCwgewogICAgY2xvbmVbZXhwZXJpbWVudCA9PSBsZW5ndGgoZXhwX2FsbG93ZWQpXQogIH0pCiAgCiAgd2ZfYnlfY3Jvc3MgPC0gd2ZfYnlfY3Jvc3NbIHdmX2J5X2Nyb3NzJGNsb25lICVpbiUgc2VsZWN0ZWRfY2xvbmVzLF0KICB3Zl9ieV9jcm9zcyA8LSAod2ZfYnlfY3Jvc3MpCiAgCiAgd2ZfYnlfY3Jvc3MkbnltcGhzIDwtIHdmX2J5X2Nyb3NzJG55bXBocysxCiAgIyBHTE0gKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqICMKICAjIEFzc3VtaW5nIE5lZ2F0aXZlIEJpbm9taWFsIGRpc3RyaWJ1dGlvbgogICMKICAKICBmaXRfbmIgPC0gTUFTUzo6Z2xtLm5iKG55bXBocyB+ICBjbG9uZSwgZGF0YSA9IHdmX2J5X2Nyb3NzKQogIEFJQyhmaXRfbmIpCgogIHBvc3Rob2MgPC0gbHNtZWFuczo6Y2xkKAogIGxzbWVhbnM6OmxzbWVhbnMoCiAgICAgZml0X25iLCB+IGNsb25lLCBhZGp1c3QgPSAidHVja2V5IiksdHlwZSA9ICJyZXNwb25zZSIKICApCiAgCiAgIyBQbG90IHBvc3QgaG9jCiAgYnJlYWtzIDwtIHBvc3Rob2MkY2xvbmVbIWdyZXBsKCJHTXxDTSIsIHBvc3Rob2MkY2xvbmUsIHBlcmwgPSBUUlVFKV0KICAKICBwbG90X2xpc3RbW3Bhc3RlKGV4cF9hbGxvd2VkKV1dIDwtIHBsb3RfY3Jvc3NfbmJfZml0KHBvc3Rob2MsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBleHBfYWxsb3dlZCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGJyZWFrcyA9IGJyZWFrcykKICAKICB9Cn0KZ2dwbG90Mjo6dGhlbWVfc2V0KGdncGxvdDI6OnRoZW1lX2dyYXkoYmFzZV9zaXplID0gMTApKQpvcmRlcmVkX2xpc3QgPC0gYyhwbG90X2xpc3RbMV0sbGlzdChlbXB0eSA9IE5VTEwpLHBsb3RfbGlzdFtjKDMsNywyLDUsNCw2KV0pCm5hbWVzKHBsb3RfbGlzdCkKbmFtZXMob3JkZXJlZF9saXN0KQpjb3dwbG90OjpwbG90X2dyaWQocGxvdGxpc3QgPSBvcmRlcmVkX2xpc3QsCiAgICAgICAgICBuY29sID0gMiwgbnJvdyA9IDQsIGFsaWduID0gInYiKQoKYGBgCg==